clear all
clc

num_cores = str2num(getenv('SLURM_CPUS_PER_TASK'));
c = parcluster;
poolobj = parpool(c, num_cores);

cd('../../../../cvx_2.2.2')
cvx_setup

warning('off','MATLAB:nargchk:deprecated')
warning('off','MATLAB:illConditionedMatrix')

cd('../bernstein_simulations/J2/codes/continuous_z')

nMC = 250;

%% DGP

J = 2;    % # of goods
N = 20000;  % # of individuals

alpha = 1;
beta = 1;

mu_x = 0;
mu_z = 0;

sigma_x = 1;
sigma_z = 1;

mu_c = -1;
sigma_c = 4;

%% PREP FOR ESTIMATION

m_regr = [4 8 4*ones(1,J-2) 8 8 4*ones(1,J-2)];

[combos_all_unique] = combinations_2(m_regr,J);

% constraints

[Aineq,bineq] = constraint_mat_theta_2(combos_all_unique,J);

Aeq = [];
beq = [];

ntheta = size(Aineq,2);

theta_NPD = zeros(ntheta,nMC);
fval_NPD = zeros(1,nMC);
flag_NPD = zeros(1,nMC);
trimmed_mean_of_ratios = zeros(nMC,1);
ratio_of_trimmed_means_second = zeros(nMC,1);
ratio_of_trimmed_means_first = zeros(nMC,1);

load('data_eval.mat')  % load data to use for evaluating ratios

parfor ss=1:nMC
    
    rng(ss)
    tic
    
    [data,data_trans] = generate_data_simult(N,J,mu_x,sigma_x,mu_z,sigma_z,alpha,beta,mu_c,sigma_c);
    
    %% ESTIMATION ROUTINE
    
    [theta_NPD(:,ss),fval_NPD(ss),flag_NPD(ss)] = estimation_NPD_cvx(data_trans,J,N,m_regr,Aineq,Aeq,bineq,beq,...
        combos_all_unique,ntheta);
        
    %% CALCULATE RATIOS

    N_temp = 10000;

    d2s_dzdz = NaN(N_temp,1);
    d2s_dzdx = NaN(N_temp,1);
    ratio_temp = NaN(N_temp,1);

    for i=1:N_temp

        [~,temp] = max(data_eval.z(i,:));

        temp_diff = setdiff(1:J,temp);

        ind_good2 = temp_diff(1);

        argument = normcdf([data_eval.x(i,temp) data_eval.x(i,temp_diff) data_eval.z(i,temp) data_eval.z(i,temp_diff)]);

        d2s_dzdz(i) = fun_d2sigma_new(theta_NPD(:,ss),argument,m_regr,combos_all_unique,J+1,J+2)*normpdf(data_eval.z(i,temp))*normpdf(data_eval.z(i,ind_good2));

        d2s_dzdx(i) = fun_d2sigma_new(theta_NPD(:,ss),argument,m_regr,combos_all_unique,J+1,2)*normpdf(data_eval.z(i,temp))*normpdf(data_eval.x(i,ind_good2));

        ratio_temp(i) = d2s_dzdz(i)/d2s_dzdx(i);

    end

    [num_temp_ord,ind_num] = sort(d2s_dzdz);
    [den_temp_ord,ind_den] = sort(d2s_dzdx);
    ratio_of_trimmed_means_second(ss) = mean(num_temp_ord(N_temp*0.25:N_temp*0.75))/mean(den_temp_ord(N_temp*0.25:N_temp*0.75));

    [ratio_temp_ord,ind] = sort(ratio_temp);
    trimmed_mean_of_ratios(ss) = mean(ratio_temp_ord(N_temp*0.25:N_temp*0.75));


    % first derivatives

    ds_dz = NaN(N_temp,1);
    ds_dx = NaN(N_temp,1);
    ratio_temp = NaN(N_temp,1);

    for i=1:N_temp

        argument = normcdf([data_eval.x(i,:) data_eval.z(i,:)]);

        ds_dz(i) = fun_dsigma(theta_NPD(:,ss),argument,m_regr,combos_all_unique,J+1)*normpdf(data_eval.z(i,1));

        ds_dx(i) = fun_dsigma(theta_NPD(:,ss),argument,m_regr,combos_all_unique,1)*normpdf(data_eval.x(i,1));

        ratio_temp(i) = ds_dz(i)/ds_dx(i);
        %i
    end

    [num_temp_ord,ind_num] = sort(ds_dz);
    [den_temp_ord,ind_den] = sort(ds_dx);
    ratio_of_trimmed_means_first(ss) = mean(num_temp_ord(N_temp*0.25:N_temp*0.75))/mean(den_temp_ord(N_temp*0.25:N_temp*0.75));

    display('*** RATIO OF TRIMMED MEANS -- SECOND ***')
    ratio_of_trimmed_means_second(ss)

    display('*** RATIO OF TRIMMED MEANS -- FIRST ***')
    ratio_of_trimmed_means_first(ss)

    ss

    toc

end

display('*** RATIO OF TRIMMED MEANS -- SECOND ***')
mean(ratio_of_trimmed_means_second)
std(ratio_of_trimmed_means_second)

display('*** RATIO OF TRIMMED MEANS -- FIRST ***')
mean(ratio_of_trimmed_means_first)
std(ratio_of_trimmed_means_first)

clear data* bineq Aineq input_list output_list

save('../../results/continuous_z/MC_simult.mat')



